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Abstract 

Starting with the Green's functions found for normal diffusion, we construct exact time- 
dependent Green's functions for subdiffusive equation with fractional time derivatives. The 
Green's functions satisfy the boundary conditions involving a linear combination of fluxes and 
concentrations. The method is particularly useful to calculate the concentration profiles in a 
multi-part system where different kinds of transport occurs in each part of it. As an example, 
we find the solution of subdiffusive equation for the system composed of two parts with normal 
diffusion and subdiffusion. 

1 Introduction 

Subdiffusion occurs in many physical systems such as amorphous semiconductors, porous media 
[1, 2] or gel solvents [3, 4]. The subdiffusion is defined by the relation 

(Ax 2 ) ~ t a (1) 

with < a < 1, where (Ace 2 ) denotes the mean square displacement of the transported particle, t 
is the time. The case of a = 1 corresponds to the normal diffusion. To obtain the concentration 
profiles (C Q ) of transported substance in subdiffusive systems one uses the subdiffusion equation 
with fractional time derivative [2, 5] 

EQ a (C a ) ee i^- t - D a ——j C a (x, t)=0, (2) 

with the Ricmann-Liouville fractional derivative defined as (for a > 0) 

d a f(t) _ 1 d n f* dt ,_ f(t') 



dt a T(n-a)dt n J (t - f )*+«-" ' 

where n is the smallest integer larger than a, D a is the subdiffusion coefficient measured in the 
units (^pr). While there is a huge literature about the normal diffusion in a variety of systems (see 
e.g. [6] and references therein), the subdiffusion is satisfactorily described in a few cases only. The 
procedure of solving the equation (2) appears to be rather complicated. Below we present a practical 
method to calculate the solutions of the equations describing in a multi-part system where diffusion 
and subdiffusion occur. The starting point is the solution corresponding to the analogous multi-part 
system where only normal diffusion is present. 
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Till now homogeneous systems have been mostly studied. However, the multi-part systems 
discussed here, in particular the membrane systems, are important for several fields of technology 
and biophysics. We mention here as examples the membranes used as micro filters and the (sub- 
diffusion processes in living cells which are crucial for the cell physiology. There are a few attempts 
described in the literature to obtain solutions of eq. (2) from the ones corresponding to the normal 
diffusion. A powerful method seems to be mapping of the solutions of diffusion equation C\(x, t) 
into the subdiffusivc ones C a (x,t) by means of the relation [7, 8, 9, 10] 

/•oo 

C a (x,t) = / n(u, t)Ci(x,u)du , (3) 
Jo 

where the Laplace transform of the function n is 

n (u, s) = 1 e -— u , (4) 

(throughout this paper the index a is assigned to the functions corresponding to subdiffusion; the 
index 1 refers to the normal diffusion). Another method is to 'rescale' the time as t — > t a [11]. Using 
the substitution one obtains the solutions of normal diffusion equation with diffusion coefficient 
depending on time. However, both procedures have been used to the homogeneous system only 
and it is unclear how to apply these approaches to a multi-part system, where the various kinds of 
diffusion occur (e.g. a porous material, where subdiffusion proceeds, is adjoined to a normal diffusive 
medium). Thus, a more general procedure is needed. 

In this paper we discuss a method of constructing solutions of the one-dimensional subdiffusive 
equation with given boundary conditions starting with those of normal diffusion, which can be 
obtained rather easily. Our method gives the Green's functions (GF) for subdiffusivc equation in 
the form of series of the Fox functions which is convenient for numerical calculations. The method 
splits into two step. In the first one the Laplace transform of the subdiffusive Green's function is 
obtained from the normal diffusion counterpart, and in the second one the Laplace transform is 
inverted. In contrast to the previous papers [7, 8, 9, 10, 11], we include boundary conditions into 
considerations. As an illustration of our method, we find the solutions of eq. (2) for the system 
consisting of two parts where the normal diffusion and subdiffusion occur, respectively. 

We restrict our considerations to the case of boundary conditions which are given in the form 
of linear combination of concentrations and fluxes. We note that this is not a serious restriction as 
the linear boundary conditions often occur in normal diffusion problems (and can be adapted for 
the subdiffusion cases). For example, for the system with fully reflecting, fully absorbing [6], and 
partially absorbing wall [12], respectively, the boundary conditions read: J{x^ t , t) = 0, C(x^ 4 ,t) = 0, 
and J(x^ M ,t) — jC(x^ I ,t), where xm denotes a position of the wall. For the partially permeable 
wall there have been used two different boundary conditions: C(xJ /[7 t) = jC(x^ [7 t) [13, 14] or 
J{x~^,t) = X(C(x^ I ,t) — C{x\ I1 t)) [15], which both should be supplemented by the continuity of 
the flux J{x~j vl ,t) = J(x^ l7 t). 



2 The method 

Since the Green's function for normal diffusion is mapped onto the subdiffusion within the method 
presented below, we first discuss the connection between the equation of normal diffusion and the 
subdiffusive equation. 



2.1 Normal diffusion 

We start with the normal diffusion equation 
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and assume that each of the boundary condition involves a linear combination of fluxes and concen- 
trations. When the system under consideration is divided into two parts by the membrane located 
at x = Xm, the boundary condition (BC) at xm is assumed to be of the form 

Li(Ci) = 6 1 C 1 (xM,t) + 6 2 C 1 (x'l I ,i) + 6 3 J 1 (xM,t) 

+ S A J 1 (xt I ,t)=0, (6) 

where the flux J\ is given as 

J 1 (x,t) = -D 1 — — — , (7) 

and Si (i = 1,2,3,4) are the parameters independent of t. The solution of the equation (5) can be 
found due to the following integral formula 



Ci{x,t)= j G 1 (x,t;x o )C 1 (xo,0)dx o , 



where Ci(x ,0) is the initial concentration, and the Green's function G\ is defined as a solution of 
the equation (5) with the initial concentration 

Gi(x, 0;x ) = 8(x - x ) , (8) 

with the appropriate boundary conditions. Since the diffusion equation is of the second order with 
respect to the space variable one needs two boundary conditions at the discontinuity of the system. 
When the system is infinite, the additional BC are needed at x — > ±oo and they are usually of the 
form G Q (x,i;xo) |a;^oo= and G a (x, t; xq) \ x ^-oo= 0, which can be treated as special cases of eq. 
(6) (8 1 — 1, 82 = 83 = 84 = for x — > —00, and 82 — 1, Si = 83 = 84 = for x — > 00, respectively). 
One usually solves the normal diffusion equation by means of the Laplace transform (LT) [6, 16] 



C(x,s) = / e' st C(x,t)dt . 
Jo 

LT of the equation (5), boundary conditions (6), and flux (7) are, respectively, 

EQ 1 (C 1 ) = sC 1 (x, s) - Di d2C *2* ,S) ~ C ^°) = . ( 9 ) 

Li(Ci) = <5i(7i(x M , s) + S 2 C 1 (x+ 1 , s) + 8 3 J 1 (x M ,s) + (5 4 Ji(x+ ,s) = , (10) 

and 

J(x,s) = -D 1 — . (11) 

The general solution of eq. (9) with the initial condition (8) is [16] 



& 1 (x,s;x ) = A 1 {s)e x ^ 7 ^ + B 1 {s)e~ x ^ 51 



I \ . p -\x-x a \y/s/Di Q2) 



where Ai(s) and Bi(s) are determined by the boundary conditions. 
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2.2 Subdiffusion 



The subdiffusive flux generated by the gradient of the concentration C a is given by the formula 

J(x t) = -D ^ dCaM 

which combined with the continuity condition 

dC a (x, t) dJ a (x, t) _ Q 
dt dx 

provides the subdiffusion equation (2). The boundary condition at x = % is assumed to be of the 
form fully analogical to that of the normal diffusion (6) i.e. 

L a (C a ) = &iC a (x M ,t) +5 2 C a (xj 4 ,t) + 8 3 J a (x M ,t) 

+ S 4 J a (x+,t)=0. (13) 

The Laplace transform of the subdiffusion equation, subdiffusive flux [17] and boundary conditions 
are, respectively, 

EQ a (C a ) = sC a (x, s) - ^d / 2 ^^ - C a (x,0) = , (14) 
J a {x,s) = -D a s^ dCa ^ s) , (15) 



L a (C a ) = 5iC a (x M ,s) + 5 2 C a (xl I ,s) + 8 3 J a (x M ,s) 
The solution of eq. (14) for the initial condition C a (x,0) = S(x — xq), is computed as 



+ <5 4 Ja(x+,s)^0. (16) 



-L \ e -|x-xo| V s°-/D a f^j\ 

where the coefficients A a (s) and B a (s) are determined, as previously, by the boundary conditions 
(16). 

2.3 From the diffusion to subdiffusion equation 

Let us note that substitution 

D x = Das 1 -* (18) 

maps the Laplace transforms of normal diffusion equation EQ\{G\) — (9), and flux J\{x, s;xo) 
(11) onto the subdiffusion ones EQ a (G a ) = (14) and J Q (x,s;xo) (15), respectively. Assuming 
that the coefficients Si, S 2 , S3, and S4 are the same for the normal diffusion and for subdiffusion 
boundary conditions (10) and (16), the substitution (18) transforms also the diffusion boundary 
condition Li(Gi) = (10) into the subdiffusion one L a (G a ) = (16). Since the Green's functions 
are uniquely determined by the subdiffusion equation, boundary conditions and initial condition 
(which are assumed to be the same for normal diffusion and subdiffusion), it is clear that the 
substitution (18) transforms the Green's function for normal diffusion (12) into the Green's function 
for subdiffusion (17). Comparing eq. (12) and eq. (17) we deduce that the coefficients A a (s) and 
B a (s) can be obtained from Ai(s) and -Bi(s) by means of the replacement (18). 
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2.4 Inverse Laplace transform 



As the Green's function is given by eq. (17), with the functions A a (s) and B a (s) of rather compli- 
cated structure (as we will see in Sec. 3.2), the calculation of the inverse Laplace transform is far 
not simple. Therefore, we expand the function into the following series 



G a (x,s;x ) = J2 a n s P+nu e x ^ sa/D « + ^ b n s<+ na e~ x ^J s * l D « 

n=0 n=0 
1 c -\x-x a \^/s<*/D a 



+ 



When the Laplace transform is in such a form, the inverse transformation can be easily calculated, 
using the formula derived by means of the method based on the Mellin transform [18] (for details of 
the derivation see the Appendix) 



wt;«)Ei- i K e -"') = 47^i (^1441 ■ ( 19 ) 



where H denotes the Fox function, a, (3 > 0, and the parameter v is not limited. Due to the formula 
(19), we obtain the Green's function as series of Fox functions. As shown in the Appendix, the 
function f v< p is given by the formula 



1 ^ sin[7r(fc / 3 + ^)]L(l + fc/3 + ^) / a\* 



sin [tt (fc/3 + v)} T (1 + k/3 + v) 

The domain of applicability of eq. (20) appears to be limited [18, 19]. The Fox function (19) 

i 

is analytical and has its series expansion (20) for every non-zero values of the variable when 
■g — 1 > 0. Otherwise, the additional restrictions should be imposed. For example, when (3 = 1, the 
equation (20) is still valid under the condition a < t. In the following we consider the cases where 
/? < 1, so we will not further discuss the domain of applicability of the formula (20). 



3 Green's functions for subdiffusive systems 

To illustrate the considerations of the previous section, we find here the Green's functions for the 
homogeneous system, which are well-known, and next we discuss a composite system containing two 
parts where, respectively, the normal diffusion and subdiffusion occur. As far as we know, such a 
system has not been considered yet. Leaving the detailed analysis of this system to the next paper, 
we derive here the respective Green's functions. 



3.1 The homogeneous system 

Since the subdiffusive homogeneous system has been already described in detail, we only show 
here that the procedure presented above gives the same results as that of the paper [2]. For the 
homogeneous system of normal diffusion, the functions A(s) and B(s) from eq. (12) vanish. So, the 
Laplace transform of GF is 

s„) = — ^e-l— °l , (21) 
and the inverse Laplace transform of (21) gives the Gausian function 
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The replacement (18) transforms the GF for normal diffusion into the following function 

G a (x,s-x ) = - * -\x-x \^JK 

and due to eq. (19) we obtain 

G ° (*•'•*•) = 2^ / «-'<»(' ;lj 7rr 1 ) 



a | x — x | \ \ £> Q i 

Let us note that the Fox function has several different representations equivalent to each other. 
Using the formula [19] 



io /V (x-x ) 2 \° , 1 A _ a i / a: - rep i f 



a\x-x \ Hl0 {\x-x \ , i-f § 



2 * 11 V VA^ 1 1 

we transform the Green's function (22) to the one derived in [2]. 
3.2 The two-part system 

As a second example, we consider a system composed of two parts: L (x < Xm) where the normal 
diffusion occurs with the diffusion coefficient D\ L and R [x > x M ) where there is the subdiffusion 
with the subdiffusion coefficient D aR . We assume the following boundary conditions at x M 

Ji(x M ,t) = J a (x M ,t) , (23) 

G 1 (x M ,t)=G a (x+,t) . (24) 

To use the procedure presented in the previous section, we temporarily assume that in the region R 
the normal diffusion occurs with the diffusion coefficient Dir. Assuming that x < Xm, the solutions 
of the diffusion equation with the boundary conditions (23) and (24) can be written in terms of the 
Laplace transform as follows 

du—(x,a;x ) = f i_ , - 1 ] e -(2x„-x-xo )x A7^T 

\2VDil8 /Dili (y/D 1L /D 1R + l) 



+ 1 p -\x-x a \^/s/D 1L (25) 



for x < xm, and 



Gu+- (X, s; X ) = 1 - e -{xM-x Q Ws/D 1L -(x-x M )^s/D 1R ^ ( 26 j 

for x > xm (the index 11 of Green's function is assigned to the system where the normal diffusion 
occurs in both parts, the index la corresponds to the system with normal diffusion and subdiffusion, 
the indices +, — refer to the signs of x — xm and xq — xm, respectively). 
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In the next step we pass from 11 to la functions. Let us note that for the homogeneous system 
the substitution (18) has been done for the whole system. However, for the composite one the 
substitution should be performed only with respect to the part where subdiffusion occurs. Thus, 

D 1R = D aR s 1 - a , (27) 

whereas the coefficient D\l remains unchanged. Substituting (27) in the functions (25) and (26) we 
obtain 

6 la —{x,8;X ) = 7 ; 1 ] e -(2x„-x-x )^7^ 

(yj D 1L s a - 1 /D aR + 1^ 

+ 1 p -\x-x \y/s/D 1L ) ( 2 g^ 

2^D 1L s 

and 

G la+ - {x, s- x ) = - — -. ,,_ e- {xM - x ^ s ' D ^-^- Xf -^ sa ' D ^ . (29) 
V D\ls + v D aR s 2 a 

It is easy to see that the functions (28) and (29) fulfil the Laplace transforms of the boundary 
conditions (23) and (24). Expanding these functions into the series, we obtain 

G la —{x,S;X ) = 1 ( e -\*-*o\y/'/DiL +e -(2x M -x-x )y/s/D 1L \ 
oo / , \ "+1 

+^=Y I -J—) s ^^ e -(^M-x-x „)V^ir 7 ( 3o) 

At i \ 1 ( I D aR \ If x M - x \ k 



n=0 fe=0 
O-gW+fc-i 



xs 



(x-XM)\fs°> /D aR 

Using the formula (19), we obtain from cqs. (30) and (31) the Green's functions for the considered 
system 

G la —(x,t;X ) = -^_( e -(-^) 2 /^ 1L t +e -(2x M -x-x Q f,AD 1L t\ 



V-DlL ^ \ V Dil J ' 2 .2^' \/Th L 



OO OC 



G la+ -(x,t;x ) = 1— 



D a R \ If X M - Xq 



n=0 k=0 

X — X M 



x /ii ^ i .^; 7 = r j. (33) 

For x > xm the calculations are very similar to the ones presented above. In this case the Green's 
functions are 



,-, / i \ 1 S^S^ D a R \ 1 f x M - x 

G la _ + (x,t;x ) = 7f -TtFT 



VDTE^ \ V d il ) k\ V V^l 
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To illustrate the above results we present two plots of the Green's functions for several values 
of a. In Fig. 1 there are presented the Green's functions (32) and (33) when the initial point Xq 
is placed at the region L (with xo = — 1). In Fig. 2 there are plots of the Green's functions (34) 
and (35) when the initial point is located at the region R (with x = 1); here Xm = 0. We use the 
arbitrary units and the values of the parameters are: D^l = 1, D\r = 0.1, D a ji = 0.1, and t = 10. 
For the double series occurring in the functions (33) and (34), and for the function f v p given by eq. 
(20), there are taken first 20 terms in each of the series. For the functions (32) and (35) given by 
the single series there are taken first 40 terms and for the function f v $ we have taken first 70 terms 
in the numerical calculations. We have checked that increasing these numbers does not change our 
results noticeably. 

Let us briefly discuss the main qualitative differences between the Green's functions for normal 
diffusion and the ones for a. < 1. For the case of x < xm the inflection point of the Green's functions 
is located at xm for a < 1 whereas for the normal diffusion this point is placed in the region R 
in a finite distance from the border between L and R. For the case of Xo > xm the differences 
between the considered Green's functions are very similar to those between the Green's function for 
a homogeneous diffusive and subdiffusive systems studied in [2]. 

4 Final remarks 

The procedure presented above allows one to solve the subdiffusion equations for the non-homogencous 
systems such as a membrane system, layered system etc. The procedure uses the fact that the normal 
diffusion equation is simpler to solve than the subdiffusive one and the relation (18) transforms the 
Green's functions, diffusive equation and boundary conditions from the normal to the subdiffusive 
case. Thus, knowing the Laplace transform of G\ one immediately obtains the Laplace transform 
ofG«. 

Summarizing our considerations, the procedure to obtain the Green's functions for the subdiffu- 
sive system with linear boundary conditions is presented as a do list: 

1. find the Laplace transforms of the Green's functions for the analogical system of the normal 
diffusion, 

2. obtain the Laplace transform of the Green's functions for subdiffusive system using the replace- 
ment (18); when the system is divided into several parts, the replacement should be applied 
only to the part where the subdiffusion occurs, 

3. expand the Laplace transform of Green's function G a into the power series of s v multiplied by 

e~ as \ 

4. use the formula (19) to calculate the inverse Laplace transform. 

The numerical calculations performed in Sec. 3.2, which give the Green's functions presented in 
Figs. 1 and 2 show usefulness of the method. It gives the explicitly analytic solutions of subdiffusive 
equations which are appropriate for numerical calculations. 
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At the end we return to the method of mapping the Green's function for normal diffusion to the 
subdiffusive one by the relation (3) which has been mentioned in the Introduction. Let us observe 
that substituting the expression (4) in cq. (3) we obtain 



-t^a 6 JO 



If we apply the replacement (18) to the r.h.s. of eq. (36), we get the identity. So, the procedure 
presented in this paper is consistent with the procedure discussed in [8] for the homogeneous system. 
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Appendix 

We derive here the relations (19) and (20). Let us consider the Laplace transform 

f^{s;a) = L[f{t)\{s) = s v e- as0 . (37) 

The Mellin transform 

M[f{x)]{p) = / xP-\f(x)dx, 
Jo 

applied to the relation (37) gives 



M 



;^( S ;a)l(p) = ^r(^). (38) 



Using the well-known formula 

M [L [/ (t)] (*)] (1 - p) = r (1 - p) M [/ (t)] (p) 

we obtain 

/ (*) = M- 1 



1 M[L[f(t)](s))(l-p) 



(39) 



r(i-p) 

Putting eq.(38) in (39) and using the inverse Mellin transform we get 

i r c+i °° i n 1 ^/^) 

2™ Jc-ioo (3a~^~ r I 1 - P) 

Comparing above equation with the definition of the Fox if-function 

/ , (oi.^i) ... (a P ,A P ) \_J_ f „ p A(p)B{p) 
np Q\ Z \ {buB,) ... (b Q ,B Q ) J-2mJ c Z C(p)D(p) dP > 

where A (p) = n£=i r - ^p), B (p) = Il"=i r (1 - <*i + (?) = ITjLi r (1 - h o + %>)> 

D (p) = Il^Li r ( a j ~ Aj'P)' an( i tne circle C separates the poles of A(p) and B(p), we obtain 



rrrnn V Res M (?) g frA 



poles ofA(p) 
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To get the poles of A(p) we take the relation bj — Bjp = —k, with k = 0, 1, 2, . . ., which gives 

bj + k 



Pj,k 

Using the relation 



lim (b j -B jP + k)T(b j -B jP )= lim J, ^ ^ + fc+1 ) 

p-p*.» v 3 JF ' y 3 JF > P ->v 3 , k {bj - Bp) {bj - B jP + 1) . . . (bj - B jP + k) 



r(i) (-i) 



k 



(-fc)(-fc + l)(-fc + 2)...(-l) " kl 

we find 

' A(p)B(p) \ _ (-l) k A' J (p^ k )B(p J . k ) 



KeSp =^{c(p)D(p)J- Bjk\C{pj, k )D( P j, k ) 



where A'j (p) = YJjT^lhpj- we obtain 

™ M _ v J_ V (- 1 ) fc - 4/ j(Pj,fc)- B (Pj,fc) . P , fc 

h ll 'h. Bjk\C{pj, k )D{ P j, k ) ' ' 

Additionally, we use the relation V (z) T (1 — z) = sin ^. z ) which gives 

1 sin (7rafc) T (1 + ak) 

T(-ak) ~ n ' 

From the above equation and the formulas (40) and (41) we get the relation (20). 
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0,18- 




Figurc 1: The Green's functions for the two-part system for xq = —1 and several values of the 
parameter a: □ corresponds to a = 0.1, otoa = 0.3, A to a = 0.5, Otoa = 0.7, *toa = 0.9, and 
V to a = 1, (the normal diffusion case). 




Figure 2: The Green's functions for the two-part system for xq = 1 and several values of the 
parameter a (the description of the symbols is the same as for Fig.l). 
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